Wellons_HW8
Madelynn Wellons
10/13/2021
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.4 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(janitor)##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testlibrary(rstan)## Loading required package: StanHeaders## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:tidyr':
##
## extractlibrary(bayesplot)## This is bayesplot version 1.8.1## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme settinglibrary(bayesrules)Exercise 6.1
a- First, you would turn the possible pi values into a finite grid (the more values, the more clear/accurate the approximation will be). Then, you would input the prior and the likelihood into the R function/into the dataset grid for this approximation in R. Then, you would calculate the likelihood and prior for each grid value, and then normalize it by dividing each value by the collective sum of values. That gives you the values for each pi “slice” that you input earlier, and you can graph this to see a visual version of the grid approximation. You can also go ahead and simulate it if you would like to check your approximation.
b- I would make more pi “slices”, as the larger the number of slices, the more accurate the approximation is. I would change it by adjusting the “values” portion of the R code in the first step to be a higher number.
Exercise 6.2 Below are my plots:
knitr::include_graphics("/Users/madelynnwellons/Pictures/HW8Plots.png")
Exercise 6.3
a- When a chain mixes slowly, the posterior will have higher levels of error, usually in increased probability in a specific range of pi values (the values that it has taken the time to “explore”)
b- When a chain has high correlation, it acts as though it is slowly mixing (it has high correlation between other values but will slowly self-correct) and so the same answer as above would apply here.
c- When a chain gets “stuck”, it is oversampling some values and so it will have a higher density at the wrong pi points (e.g. higher density than there should be on the tail of a distribution)
Exercise 6.4
a- It is important to look at MCMC diagnostics to ensure that there were no issues with the chains, as if there were then the posterior model would be very off.
b- MCMC simulations are helpful because they can approximate posteriors quickly that are impossible to calculate by hand and would take hours/days to simulate properly.
c- RStan allows you to use Stan code in order to create the data model that the MCMC will be using.
d- I understand most of the chapter, I am just having trouble figuring out when we would actually use this- Nico mentioned that we wouldn’t since there are better approximation methods out there now, but why did people use this in the first place?
Exercise 6.5
a- First we will check what the posterior should be from an analytical standpoint:
summarize_beta_binomial(3, 8, 2, 10)## model alpha beta mean mode var sd
## 1 prior 3 8 0.2727273 0.2222222 0.016528926 0.12856487
## 2 posterior 5 16 0.2380952 0.2105263 0.008245723 0.09080596We know that the posterior should be (5, 16), this will come in handy later.
Now we will split up the pi values in the grid:
grid_data65 <- data.frame(pi_grid = seq(from = 0, to = 1, length = 5))Then we will use dbeta and dbinom to evaluate the prior and likelihood given the data:
grid_data65 <- grid_data65 |>
mutate(prior = dbeta(pi_grid, 3, 8),
likelihood = dbinom(2, 10, pi_grid))#Step 3: approximate the posterior
grid_data65 <- grid_data65 |>
mutate(unnormalized = likelihood * prior,
posterior = unnormalized / sum(unnormalized))
#Confirming that posterior approximation sums to 1
grid_data65 |>
summarize(sum(unnormalized), sum(posterior))## sum(unnormalized) sum(posterior)
## 1 0.8765603 1#Examine the grid approx. posterior
round(grid_data65, 2)## pi_grid prior likelihood unnormalized posterior
## 1 0.00 0.00 0.00 0.00 0.00
## 2 0.25 3.00 0.28 0.85 0.96
## 3 0.50 0.70 0.04 0.03 0.04
## 4 0.75 0.01 0.00 0.00 0.00
## 5 1.00 0.00 0.00 0.00 0.00#Plot grid approx. posterior
ggplot(grid_data65, aes(x = pi_grid, y = posterior)) +
geom_point() +
geom_segment(aes(x = pi_grid, xend = pi_grid, y = 0, yend = posterior))
#Sample from discretized posterior
set.seed(3456)
post_sample65 <- sample_n(grid_data65, size = 10000,
weight = posterior, replace = TRUE)
post_sample65 |>
tabyl(pi_grid) |>
adorn_totals("row")## pi_grid n percent
## 0.25 9651 0.9651
## 0.5 348 0.0348
## 0.75 1 0.0001
## Total 10000 1.0000#Histogram of grid simulation with posterior pdf
ggplot(post_sample65, aes(x = pi_grid)) +
geom_histogram(aes(y = ..density..), color = "white") +
stat_function(fun = dbeta, args = list(5, 16)) +
lims(x = c(0, 1))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing missing values (geom_bar).
b- We can likely better approximate the posterior with more “slices”- in the graph above, you can see how only having the points for a few select pi values makes it not as accurate when compared to the real posterior.
Now we will split up the pi values in the grid:
grid_data65b <- data.frame(pi_grid = seq(from = 0, to = 1, length = 201))Then we will use dbeta and dbinom to evaluate the prior and likelihood given the data:
grid_data65b <- grid_data65b |>
mutate(prior = dbeta(pi_grid, 3, 8),
likelihood = dbinom(2, 10, pi_grid))#Step 3: approximate the posterior
grid_data65b <- grid_data65b |>
mutate(unnormalized = likelihood * prior,
posterior = unnormalized / sum(unnormalized))
#Confirming that posterior approximation sums to 1
grid_data65b |>
summarize(sum(unnormalized), sum(posterior))## sum(unnormalized) sum(posterior)
## 1 41.79567 1#Plot grid approx. posterior
ggplot(grid_data65b, aes(x = pi_grid, y = posterior)) +
geom_point() +
geom_segment(aes(x = pi_grid, xend = pi_grid, y = 0, yend = posterior))
set.seed(3456)
post_sample65b <- sample_n(grid_data65b, size = 10000,
weight = posterior, replace = TRUE)
post_sample65b |>
tabyl(pi_grid) |>
adorn_totals("row")## pi_grid n percent
## 0.02 1 0.0001
## 0.03 3 0.0003
## 0.035 7 0.0007
## 0.04 5 0.0005
## 0.045 11 0.0011
## 0.05 7 0.0007
## 0.055 9 0.0009
## 0.06 15 0.0015
## 0.065 21 0.0021
## 0.07 40 0.0040
## 0.075 37 0.0037
## 0.08 43 0.0043
## 0.085 55 0.0055
## 0.09 58 0.0058
## 0.095 80 0.0080
## 0.1 80 0.0080
## 0.105 79 0.0079
## 0.11 109 0.0109
## 0.115 117 0.0117
## 0.12 112 0.0112
## 0.125 135 0.0135
## 0.13 131 0.0131
## 0.135 127 0.0127
## 0.14 153 0.0153
## 0.145 186 0.0186
## 0.15 175 0.0175
## 0.155 200 0.0200
## 0.16 194 0.0194
## 0.165 191 0.0191
## 0.17 203 0.0203
## 0.175 187 0.0187
## 0.18 212 0.0212
## 0.185 183 0.0183
## 0.19 174 0.0174
## 0.195 191 0.0191
## 0.2 211 0.0211
## 0.205 225 0.0225
## 0.21 230 0.0230
## 0.215 212 0.0212
## 0.22 219 0.0219
## 0.225 205 0.0205
## 0.23 193 0.0193
## 0.235 214 0.0214
## 0.24 210 0.0210
## 0.245 212 0.0212
## 0.25 189 0.0189
## 0.255 220 0.0220
## 0.26 183 0.0183
## 0.265 183 0.0183
## 0.27 193 0.0193
## 0.275 198 0.0198
## 0.28 183 0.0183
## 0.285 193 0.0193
## 0.29 173 0.0173
## 0.295 135 0.0135
## 0.3 150 0.0150
## 0.305 136 0.0136
## 0.31 128 0.0128
## 0.315 142 0.0142
## 0.32 112 0.0112
## 0.325 118 0.0118
## 0.33 116 0.0116
## 0.335 101 0.0101
## 0.34 117 0.0117
## 0.345 84 0.0084
## 0.35 95 0.0095
## 0.355 97 0.0097
## 0.36 75 0.0075
## 0.365 100 0.0100
## 0.37 62 0.0062
## 0.375 78 0.0078
## 0.38 60 0.0060
## 0.385 61 0.0061
## 0.39 64 0.0064
## 0.395 46 0.0046
## 0.4 28 0.0028
## 0.405 50 0.0050
## 0.41 36 0.0036
## 0.415 45 0.0045
## 0.42 36 0.0036
## 0.425 40 0.0040
## 0.43 32 0.0032
## 0.435 23 0.0023
## 0.44 24 0.0024
## 0.445 27 0.0027
## 0.45 14 0.0014
## 0.455 19 0.0019
## 0.46 19 0.0019
## 0.465 21 0.0021
## 0.47 17 0.0017
## 0.475 9 0.0009
## 0.48 15 0.0015
## 0.485 7 0.0007
## 0.49 7 0.0007
## 0.495 10 0.0010
## 0.5 10 0.0010
## 0.505 5 0.0005
## 0.51 9 0.0009
## 0.515 6 0.0006
## 0.52 3 0.0003
## 0.525 3 0.0003
## 0.53 6 0.0006
## 0.535 5 0.0005
## 0.54 2 0.0002
## 0.545 3 0.0003
## 0.55 1 0.0001
## 0.555 3 0.0003
## 0.57 1 0.0001
## 0.575 1 0.0001
## 0.58 1 0.0001
## 0.585 1 0.0001
## 0.59 1 0.0001
## 0.6 1 0.0001
## 0.605 1 0.0001
## 0.61 1 0.0001
## 0.62 1 0.0001
## 0.63 1 0.0001
## 0.735 1 0.0001
## Total 10000 1.0000#Histogram of grid simulation with posterior pdf
ggplot(post_sample65b, aes(x = pi_grid)) +
geom_histogram(aes(y = ..density..), color = "white") +
stat_function(fun = dbeta, args = list(5, 16)) +
lims(x = c(0, 1))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing missing values (geom_bar).
Exercise 6.6
a- We can follow similar steps as we did in the last problem, but adjusting it for Gamma-Poisson.
# Step 1: Define a grid
grid_data6 <- data.frame(lambda_grid = seq(from = 0, to = 8, length = 9))
# Step 2: Evaluate the prior & likelihood at each lambda
grid_data6 <- grid_data6 |>
mutate(prior = dgamma(lambda_grid, 20, 5),
likelihood = dpois(0, lambda_grid) * dpois(1, lambda_grid) * dpois(0, lambda_grid))
# Step 3: Approximate the posterior
grid_data6 <- grid_data6 |>
mutate(unnormalized = likelihood * prior,
posterior = unnormalized / sum(unnormalized))
# Set the seed
set.seed(4635)
# Step 4: sample from the discretized posterior
post_sample6 <- sample_n(grid_data6, size = 10000,
weight = posterior, replace = TRUE)Now we can check what the value of the posterior should be analytically:
summarize_gamma_poisson(20, 5, 1, 3)## model shape rate mean mode var sd
## 1 prior 20 5 4.000 3.8 0.800000 0.8944272
## 2 posterior 21 8 2.625 2.5 0.328125 0.5728220Now we can plot this on the histogram:
# Histogram of the grid simulation with posterior pdf
ggplot(post_sample6, aes(x = lambda_grid)) +
geom_histogram(aes(y = ..density..), color = "white") +
stat_function(fun = dgamma, args = list(21, 8)) +
lims(x = c(0, 8))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing missing values (geom_bar).
b- Now we can try that all again but with 201 grid values!
# Step 1: Define a grid
grid_data6b <- data.frame(lambda_grid = seq(from = 0, to = 8, length = 201))
# Step 2: Evaluate the prior & likelihood at each lambda
grid_data6b <- grid_data6b |>
mutate(prior = dgamma(lambda_grid, 20, 5),
likelihood = dpois(0, lambda_grid) * dpois(1, lambda_grid) * dpois(0, lambda_grid))
# Step 3: Approximate the posterior
grid_data6b <- grid_data6b |>
mutate(unnormalized = likelihood * prior,
posterior = unnormalized / sum(unnormalized))
# Set the seed
set.seed(4635)
# Step 4: sample from the discretized posterior
post_sample6b <- sample_n(grid_data6b, size = 10000,
weight = posterior, replace = TRUE)Now we can plot this on the histogram:
# Histogram of the grid simulation with posterior pdf
ggplot(post_sample6b, aes(x = lambda_grid)) +
geom_histogram(aes(y = ..density..), color = "white") +
stat_function(fun = dgamma, args = list(21, 8)) +
lims(x = c(0, 8))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing missing values (geom_bar).
Exercise 6.7
a- Here we will once again use the same methods to create the grid approximation (but adjusting the code for it being a normal-normal):
# Step 1: Define a grid
grid_data7 <- data.frame(mu_grid = seq(from = 5, to = 15, length = 11))# Step 2: Evaluate the prior & likelihood at each mu- similar to last two but code differs as the mean is the first value in the prior and the sd is the second, and similarly for the posterior
grid_data7 <- grid_data7 |>
mutate(prior = dnorm(mu_grid, mean = 10, sd = 1.2),
likelihood = dnorm(7.1, mean = mu_grid, sd = 1.3)*
dnorm(8.9, mean = mu_grid, sd = 1.3)*
dnorm(8.4, mean = mu_grid, sd = 1.3)*
dnorm(8.6, mean = mu_grid, sd = 1.3))#Step 3: approximate the posterior, this is almost identical to our usual code here
grid_data7 <- grid_data7 %>%
mutate(unnormalized = likelihood * prior,
posterior = unnormalized / sum(unnormalized))
ggplot(grid_data7, aes(x = mu_grid, y = posterior)) +
geom_point() +
geom_segment(aes(x = mu_grid, xend = mu_grid, y = 0, yend = posterior))
I now realize that we didn’t actually have to do the simulations for the last few questions as the code that was given to us for the normal-normal didn’t include that, sorry that there is extra code to review!!
b- Now we can do the same thing, but with 201 mu values!
# Step 1: Define a grid
grid_data7b <- data.frame(mu_grid = seq(from = 5, to = 15, length = 201))
# Step 2: Evaluate the prior & likelihood at each mu
grid_data7b <- grid_data7b |>
mutate(prior = dnorm(mu_grid, mean = 10, sd = 1.2),
likelihood = dnorm(7.1, mean = mu_grid, sd = 1.3)*
dnorm(8.9, mean = mu_grid, sd = 1.3)*
dnorm(8.4, mean = mu_grid, sd = 1.3)*
dnorm(8.6, mean = mu_grid, sd = 1.3))
#Step 3: approximate the posterior
grid_data7b <- grid_data7b %>%
mutate(unnormalized = likelihood * prior,
posterior = unnormalized / sum(unnormalized))
ggplot(grid_data7b, aes(x = mu_grid, y = posterior)) +
geom_point() +
geom_segment(aes(x = mu_grid, xend = mu_grid, y = 0, yend = posterior))
Exercise 6.8
a- One situation in which we might want inference for multiple parameters is if you are wanting to estimate both the height and weight of a child based on their age (height being one variable, weight being another).
b- Dimensionality adds another layer (“dimension”) or axis to the grid approximation; this causes us to need an incredibly large amount of pi/lambda/mu “slices” in order to get a decent approximation of the posterior (e.g. if you need 100 on each axis, you’d now need to calculate 10000 values in this approximation, whereas if it was one-dimensional you would only need 100). This is a curse because it makes it nearly impossible to get a good grid approximation without the analysis taking too long.
Exercise 6.9
a- Both MCMC and grid approximation are approximations of the true posterior (not an exact analytical calculation)
b- Both MCMC and grid approximation can approximate the posterior for those that we cannot calculate analytically (or at least it would be very difficult to do so)
c- Grid approximation has the advantage of each value being independent from another, and that reduces the likelihood of error that comes with MCMC
d- MCMC has the advantage of being able to calculate more complex posteriors in a very quick amount of time compared to grid approximation.
Exercise 6.10
a- Yes, this is a Markov chain; each day is dependent on the last (e.g. you would likely not order it two days in a row)
b- No, this is not a Markov chain, as the values are independent
c- Yes, this is a Markov chain; each day you learn more from your chess matches with your roommate (and vice versa) so each day is dependent on the rest.
Exercise 6.11
a- Below is the code for the rstan model for this problem:
bb_model <- "
data {
int<lower = 0, upper = 20> Y;
}
parameters {
real<lower = 0, upper = 1> pi;
}
model {
Y ~ binomial(20, pi);
pi ~ beta(1, 1);
}
"b-
#I have no value here for the number of observations since it was not given in the text
gp_model <- "
data {
int<lower = 0> Y;
}
parameters {
real<lower = 0> lambda;
}
model {
Y ~ poisson(lambda);
lambda ~ gamma(4, 2);
}
"c-
normal_model <- '
data {
vector[] Y; //Once again the value was not given so I am leaving the vector[] blank here
}
parameters {
real mu;
}
model {
Y ~ normal(mu, 1); //This is the (mu, 1^2) we were given in the problem
mu ~ normal(0, 10); //This is the (0, 10^2) prior we were given in the problem
}
'Exercise 6.12
a-
#Step 1: defining the model
bb_model <- "
data {
int<lower = 0, upper = 20> Y;
}
parameters {
real<lower = 0, upper = 1> pi;
}
model {
Y ~ binomial(20, pi);
pi ~ beta(1, 1);
}
"
#Step 2: simulate the posterior
bb_sim <- stan(model_code = bb_model, data = list(Y = 12),
chains = 4, iter = 5000*2, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL '79757aaadd38925ba0ee929aa4f540cd' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.035405 seconds (Warm-up)
## Chain 1: 0.037459 seconds (Sampling)
## Chain 1: 0.072864 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '79757aaadd38925ba0ee929aa4f540cd' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.035106 seconds (Warm-up)
## Chain 2: 0.03768 seconds (Sampling)
## Chain 2: 0.072786 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '79757aaadd38925ba0ee929aa4f540cd' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.036579 seconds (Warm-up)
## Chain 3: 0.038127 seconds (Sampling)
## Chain 3: 0.074706 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '79757aaadd38925ba0ee929aa4f540cd' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.036117 seconds (Warm-up)
## Chain 4: 0.038943 seconds (Sampling)
## Chain 4: 0.07506 seconds (Total)
## Chain 4:b- The model here will be the same as in 6.11, but we will have defined the number of observations since we only have one datapoint for Y:
#Step 1: define the model
gp_model <- "
data {
int<lower = 0> Y[1];
}
parameters {
real<lower = 0> lambda;
}
model {
Y ~ poisson(lambda);
lambda ~ gamma(4, 2);
}
"
#Step 2: simulate the posterior
gp_sim <- stan(model_code = gp_model, data = list(Y = c(3)),
chains = 4, iter = 5000*2, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Error in mod$fit_ptr() :
## Exception: mismatch in number dimensions declared and found in context; processing stage=data initialization; variable name=Y; dims declared=(1); dims found=() (in 'model744c75b14bad_5d7114e32bd0158afb619c8531758230' at line 3)## failed to create the sampler; sampling not donec- The model here will be the same as in 6.11, but we will have defined the number of observations since we only have one datapoint for Y:
#step 1: define the model
normal_model <- '
data {
vector[1] Y;
}
parameters {
real mu;
}
model {
Y ~ normal(mu, 1); //This is the (mu, 1^2) we were given in the problem
mu ~ normal(0, 10); //This is the (0, 10^2) prior we were given in the problem
}
'
#Step 2: simulate the posterior
gn_sim <- stan(model_code = normal_model, data = list(Y = c(12.2)), #Here I adjusted Nico's code since he made this based on the earlier problem with 4 values- I left the data in since it was only one value and did not create a vector like he did
chains = 4, iter = 5000*2, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Error in mod$fit_ptr() :
## Exception: mismatch in number dimensions declared and found in context; processing stage=data initialization; variable name=Y; dims declared=(1); dims found=() (in 'model744c208852b1_4deb2deaf6c6dd8a4d87f15451909dea' at line 3)## failed to create the sampler; sampling not doneExercise 6.13
a-
#Step 1: defining the model
bb_model <- "
data {
int<lower = 0, upper = 10> Y;
}
parameters {
real<lower = 0, upper = 1> pi;
}
model {
Y ~ binomial(10, pi);
pi ~ beta(3, 8);
}
"
#Step 2: simulate the posterior
bb_sim <- stan(model_code = bb_model, data = list(Y = 2),
chains = 3, iter = 12000, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL '6cc1e5f794ef111dd783f9c98fed9c41' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.7e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 12000 [ 0%] (Warmup)
## Chain 1: Iteration: 1200 / 12000 [ 10%] (Warmup)
## Chain 1: Iteration: 2400 / 12000 [ 20%] (Warmup)
## Chain 1: Iteration: 3600 / 12000 [ 30%] (Warmup)
## Chain 1: Iteration: 4800 / 12000 [ 40%] (Warmup)
## Chain 1: Iteration: 6000 / 12000 [ 50%] (Warmup)
## Chain 1: Iteration: 6001 / 12000 [ 50%] (Sampling)
## Chain 1: Iteration: 7200 / 12000 [ 60%] (Sampling)
## Chain 1: Iteration: 8400 / 12000 [ 70%] (Sampling)
## Chain 1: Iteration: 9600 / 12000 [ 80%] (Sampling)
## Chain 1: Iteration: 10800 / 12000 [ 90%] (Sampling)
## Chain 1: Iteration: 12000 / 12000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.041625 seconds (Warm-up)
## Chain 1: 0.040675 seconds (Sampling)
## Chain 1: 0.0823 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '6cc1e5f794ef111dd783f9c98fed9c41' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 12000 [ 0%] (Warmup)
## Chain 2: Iteration: 1200 / 12000 [ 10%] (Warmup)
## Chain 2: Iteration: 2400 / 12000 [ 20%] (Warmup)
## Chain 2: Iteration: 3600 / 12000 [ 30%] (Warmup)
## Chain 2: Iteration: 4800 / 12000 [ 40%] (Warmup)
## Chain 2: Iteration: 6000 / 12000 [ 50%] (Warmup)
## Chain 2: Iteration: 6001 / 12000 [ 50%] (Sampling)
## Chain 2: Iteration: 7200 / 12000 [ 60%] (Sampling)
## Chain 2: Iteration: 8400 / 12000 [ 70%] (Sampling)
## Chain 2: Iteration: 9600 / 12000 [ 80%] (Sampling)
## Chain 2: Iteration: 10800 / 12000 [ 90%] (Sampling)
## Chain 2: Iteration: 12000 / 12000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.040583 seconds (Warm-up)
## Chain 2: 0.051114 seconds (Sampling)
## Chain 2: 0.091697 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '6cc1e5f794ef111dd783f9c98fed9c41' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 8e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 12000 [ 0%] (Warmup)
## Chain 3: Iteration: 1200 / 12000 [ 10%] (Warmup)
## Chain 3: Iteration: 2400 / 12000 [ 20%] (Warmup)
## Chain 3: Iteration: 3600 / 12000 [ 30%] (Warmup)
## Chain 3: Iteration: 4800 / 12000 [ 40%] (Warmup)
## Chain 3: Iteration: 6000 / 12000 [ 50%] (Warmup)
## Chain 3: Iteration: 6001 / 12000 [ 50%] (Sampling)
## Chain 3: Iteration: 7200 / 12000 [ 60%] (Sampling)
## Chain 3: Iteration: 8400 / 12000 [ 70%] (Sampling)
## Chain 3: Iteration: 9600 / 12000 [ 80%] (Sampling)
## Chain 3: Iteration: 10800 / 12000 [ 90%] (Sampling)
## Chain 3: Iteration: 12000 / 12000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.040984 seconds (Warm-up)
## Chain 3: 0.045898 seconds (Sampling)
## Chain 3: 0.086882 seconds (Total)
## Chain 3:b-
#Creating trace plot
mcmc_trace(bb_sim, pars = "pi", size = .1)
c- The x-axis has values of 0 to 6000; it only goes to 6000 instead of 12000 because it removed the first half of the values as the “burn in” period.
d- Below are the density plots for individual chains:
mcmc_dens_overlay(bb_sim, pars = "pi") +
ylab("density")
e- Below is the calculation of the posterior model:
#alpha posterior = alpha prior + y
3+2## [1] 5#beta posterior = beta prior + n - y
8+10-2## [1] 16The posterior here should be (5, 16); I am unsure of how to plot the MCMC density plot over the posterior, but I will plot the posterior below so we can compare:
plot_beta_binomial(3, 8, 2, 10)
The posterior here does look similar to our density plot!
Exercise 6.14
a-
#Step 1: defining the model
bb_model <- "
data {
int<lower = 0, upper = 12> Y;
}
parameters {
real<lower = 0, upper = 1> pi;
}
model {
Y ~ binomial(12, pi);
pi ~ beta(4, 3);
}
"
#Step 2: simulate the posterior
bb_sim <- stan(model_code = bb_model, data = list(Y = 4),
chains = 3, iter = 12000, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL '981ca33f2a8a2fb571d6627286410841' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.32 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 12000 [ 0%] (Warmup)
## Chain 1: Iteration: 1200 / 12000 [ 10%] (Warmup)
## Chain 1: Iteration: 2400 / 12000 [ 20%] (Warmup)
## Chain 1: Iteration: 3600 / 12000 [ 30%] (Warmup)
## Chain 1: Iteration: 4800 / 12000 [ 40%] (Warmup)
## Chain 1: Iteration: 6000 / 12000 [ 50%] (Warmup)
## Chain 1: Iteration: 6001 / 12000 [ 50%] (Sampling)
## Chain 1: Iteration: 7200 / 12000 [ 60%] (Sampling)
## Chain 1: Iteration: 8400 / 12000 [ 70%] (Sampling)
## Chain 1: Iteration: 9600 / 12000 [ 80%] (Sampling)
## Chain 1: Iteration: 10800 / 12000 [ 90%] (Sampling)
## Chain 1: Iteration: 12000 / 12000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.043275 seconds (Warm-up)
## Chain 1: 0.0454 seconds (Sampling)
## Chain 1: 0.088675 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '981ca33f2a8a2fb571d6627286410841' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 12000 [ 0%] (Warmup)
## Chain 2: Iteration: 1200 / 12000 [ 10%] (Warmup)
## Chain 2: Iteration: 2400 / 12000 [ 20%] (Warmup)
## Chain 2: Iteration: 3600 / 12000 [ 30%] (Warmup)
## Chain 2: Iteration: 4800 / 12000 [ 40%] (Warmup)
## Chain 2: Iteration: 6000 / 12000 [ 50%] (Warmup)
## Chain 2: Iteration: 6001 / 12000 [ 50%] (Sampling)
## Chain 2: Iteration: 7200 / 12000 [ 60%] (Sampling)
## Chain 2: Iteration: 8400 / 12000 [ 70%] (Sampling)
## Chain 2: Iteration: 9600 / 12000 [ 80%] (Sampling)
## Chain 2: Iteration: 10800 / 12000 [ 90%] (Sampling)
## Chain 2: Iteration: 12000 / 12000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.040989 seconds (Warm-up)
## Chain 2: 0.047734 seconds (Sampling)
## Chain 2: 0.088723 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '981ca33f2a8a2fb571d6627286410841' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 12000 [ 0%] (Warmup)
## Chain 3: Iteration: 1200 / 12000 [ 10%] (Warmup)
## Chain 3: Iteration: 2400 / 12000 [ 20%] (Warmup)
## Chain 3: Iteration: 3600 / 12000 [ 30%] (Warmup)
## Chain 3: Iteration: 4800 / 12000 [ 40%] (Warmup)
## Chain 3: Iteration: 6000 / 12000 [ 50%] (Warmup)
## Chain 3: Iteration: 6001 / 12000 [ 50%] (Sampling)
## Chain 3: Iteration: 7200 / 12000 [ 60%] (Sampling)
## Chain 3: Iteration: 8400 / 12000 [ 70%] (Sampling)
## Chain 3: Iteration: 9600 / 12000 [ 80%] (Sampling)
## Chain 3: Iteration: 10800 / 12000 [ 90%] (Sampling)
## Chain 3: Iteration: 12000 / 12000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.043588 seconds (Warm-up)
## Chain 3: 0.047709 seconds (Sampling)
## Chain 3: 0.091297 seconds (Total)
## Chain 3:b-
#Creating trace plot
mcmc_trace(bb_sim, pars = "pi", size = .1)
c- The x-axis has values of 0 to 6000; it only goes to 6000 instead of 12000 because it removed the first half of the values as the “burn in” period.
d- Below are the density plots for individual chains:
mcmc_dens_overlay(bb_sim, pars = "pi") +
ylab("density")
#alpha posterior = alpha prior + y
4+4## [1] 8#beta posterior = beta prior + n - y
3+12-4## [1] 11The posterior here should be (8, 11); once again, I am unsure of how to plot the MCMC density plot over the posterior, but I will plot the posterior below so we can compare:
plot_beta_binomial(4, 3, 4, 12)
This also looks quite similar, although it looks like Chain 1 was the closest (as it had a slightly higher peak than the other 2 chains).
Exercise 6.15
a-
#Step 1: define the model
gp_model <- "
data {
int<lower = 0> Y[3];
}
parameters {
real<lower = 0> lambda;
}
model {
Y ~ poisson(lambda);
lambda ~ gamma(20, 5);
}
"
#Step 2: simulate the posterior
gp_sim <- stan(model_code = gp_model, data = list(Y = c(0, 1, 0)),
chains = 4, iter = 10000, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'f65cbc91cb99a701a7ac1a425b95cf50' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.031433 seconds (Warm-up)
## Chain 1: 0.035804 seconds (Sampling)
## Chain 1: 0.067237 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'f65cbc91cb99a701a7ac1a425b95cf50' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 8e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.032075 seconds (Warm-up)
## Chain 2: 0.032144 seconds (Sampling)
## Chain 2: 0.064219 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'f65cbc91cb99a701a7ac1a425b95cf50' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.032049 seconds (Warm-up)
## Chain 3: 0.031272 seconds (Sampling)
## Chain 3: 0.063321 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'f65cbc91cb99a701a7ac1a425b95cf50' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 7e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.030739 seconds (Warm-up)
## Chain 4: 0.031391 seconds (Sampling)
## Chain 4: 0.06213 seconds (Total)
## Chain 4:b-
#Creating trace plot
mcmc_trace(gp_sim, pars = "lambda", size = .1)
#density plot for the 4 chains
mcmc_dens_overlay(gp_sim, pars = "lambda") +
ylab("density")
c- The most posterior plausible value of lambda appears to be about 2.6 or 2.7
d- We can find the exact posterior using the summarize_gamma_poisson function:
summarize_gamma_poisson(20, 5, 1, 3)## model shape rate mean mode var sd
## 1 prior 20 5 4.000 3.8 0.800000 0.8944272
## 2 posterior 21 8 2.625 2.5 0.328125 0.5728220My approximation seems fairly close! The estimated mean and mode are close to my estimate of the most likely value of lambda.
Exercise 6.16
a-
#Step 1: define the model
gp_model <- "
data {
int<lower = 0> Y[3];
}
parameters {
real<lower = 0> lambda;
}
model {
Y ~ poisson(lambda);
lambda ~ gamma(5, 5);
}
"
#Step 2: simulate the posterior
gp_sim <- stan(model_code = gp_model, data = list(Y = c(0, 1, 0)),
chains = 4, iter = 10000, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'd6d79c0ae95ef30a3463fe7b9e36388e' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.032599 seconds (Warm-up)
## Chain 1: 0.034813 seconds (Sampling)
## Chain 1: 0.067412 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'd6d79c0ae95ef30a3463fe7b9e36388e' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.033161 seconds (Warm-up)
## Chain 2: 0.03648 seconds (Sampling)
## Chain 2: 0.069641 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'd6d79c0ae95ef30a3463fe7b9e36388e' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.032811 seconds (Warm-up)
## Chain 3: 0.033501 seconds (Sampling)
## Chain 3: 0.066312 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'd6d79c0ae95ef30a3463fe7b9e36388e' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 7e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.034918 seconds (Warm-up)
## Chain 4: 0.032646 seconds (Sampling)
## Chain 4: 0.067564 seconds (Total)
## Chain 4:b-
#Creating trace plot
mcmc_trace(gp_sim, pars = "lambda", size = .1)
#density plot for the 4 chains
mcmc_dens_overlay(gp_sim, pars = "lambda") +
ylab("density")
c- The most plausible posterior value of lambda seems to be about .6!
d- We can once again use the summarize function to determine the posterior model:
summarize_gamma_poisson(5, 5, 1, 3)## model shape rate mean mode var sd
## 1 prior 5 5 1.00 0.800 0.20000 0.4472136
## 2 posterior 6 8 0.75 0.625 0.09375 0.3061862Here the posterior is (6, 8) and it appears to be relatively close to my MCMC approximation based on the mean and mode values!
Exercise 6.17
a-
normal_model <- '
data {
vector[4] Y;
}
parameters {
real mu;
}
model {
Y ~ normal(mu, 1.3);
mu ~ normal(10, 1.2);
}
'
d <- list(Y = c(7.1, 8.9, 8.4, 8.6)) #this time since there is more than one observation I will do what Nico did and put the values into a vector
gn_sim <- stan(model_code = normal_model, data = d,
chains = 4, iter = 10000, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'aafc09a58798752b5f677b912ca25e71' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.031024 seconds (Warm-up)
## Chain 1: 0.036436 seconds (Sampling)
## Chain 1: 0.06746 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'aafc09a58798752b5f677b912ca25e71' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.030576 seconds (Warm-up)
## Chain 2: 0.034138 seconds (Sampling)
## Chain 2: 0.064714 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'aafc09a58798752b5f677b912ca25e71' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.032101 seconds (Warm-up)
## Chain 3: 0.031564 seconds (Sampling)
## Chain 3: 0.063665 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'aafc09a58798752b5f677b912ca25e71' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.030041 seconds (Warm-up)
## Chain 4: 0.033992 seconds (Sampling)
## Chain 4: 0.064033 seconds (Total)
## Chain 4:b-
#Creating trace plot
mcmc_trace(gn_sim, pars = "mu", size = .1)
#density plot for the 4 chains
mcmc_dens_overlay(gn_sim, pars = "mu") +
ylab("density")
c- The most psterior plausible value of mu appears to be about 8.7
d- Here we can use the summarize normal normal function to calculate the posterior:
summarize_normal_normal(10, 1.2, 1.3, mean(c(7.1, 8.9, 8.4, 8.6)), 4)## model mean mode var sd
## 1 prior 10.00000 10.00000 1.4400000 1.2000000
## 2 posterior 8.64698 8.64698 0.3266577 0.5715398Once again, this also seems very close to my MCMC approximation based on the mean/mode! Here the distribution would be (8.65, .57) rounded to 2 decimal points.
Exercise 6.18
a-
normal_model <- '
data {
vector[5] Y;
}
parameters {
real mu;
}
model {
Y ~ normal(mu, 8);
mu ~ normal(-14, 2);
}
'
d <- list(Y = c(-10.1, 5.5, .1, -1.4, 11.5)) #this time since there is more than one observation I will do what Nico did and put the values into a vector
gn_sim <- stan(model_code = normal_model, data = d,
chains = 4, iter = 10000, seed = 1234)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'cdfc0c053324bc0fce964d8c0a25d259' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.032517 seconds (Warm-up)
## Chain 1: 0.034787 seconds (Sampling)
## Chain 1: 0.067304 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'cdfc0c053324bc0fce964d8c0a25d259' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.031237 seconds (Warm-up)
## Chain 2: 0.032545 seconds (Sampling)
## Chain 2: 0.063782 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'cdfc0c053324bc0fce964d8c0a25d259' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.030989 seconds (Warm-up)
## Chain 3: 0.032268 seconds (Sampling)
## Chain 3: 0.063257 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'cdfc0c053324bc0fce964d8c0a25d259' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 7e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.031376 seconds (Warm-up)
## Chain 4: 0.034067 seconds (Sampling)
## Chain 4: 0.065443 seconds (Total)
## Chain 4:b-
#Creating trace plot
mcmc_trace(gn_sim, pars = "mu", size = .1)
#density plot for the 4 chains
mcmc_dens_overlay(gn_sim, pars = "mu") +
ylab("density")
c- Here the most likely posterior value of mu appears to be roughly -11 or -10.
d- Below is the calculated posterior:
summarize_normal_normal(-14, 2, 8, mean(c(-10.1, 5.5, .1, -1.4, 11.5)), 5)## model mean mode var sd
## 1 prior -14.0 -14.0 4.000000 2.000000
## 2 posterior -10.4 -10.4 3.047619 1.745743The posterior here is (-10.4, 1.75) rounded to 2 decimal points, and appears to be very close to our approximated posterior based on the mean/mode!